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Abstract 

The spectral properties of Kitaev's honeycomb lattice model are investigated both 
analytically and numerically with the focus on the non-abelian phase of the model. 
After summarizing the fermionization technique which maps spins into free Majo- 
rana fermions, we evaluate the spectrum of sparse vortex configurations and derive 
the interaction between two vortices as a function of their separation. We consider 
the effect vortices can have on the fermionic spectrum as well as on the phase tran- 
sition between the abelian and non-abelian phases. We explicitly demonstrate the 
2"-fold ground state degeneracy in the presence of 2n well separated vortices and 
the lifting of the degeneracy due to their short-range interactions. The calculations 
are performed on an infinite lattice. In addition to the analytic treatment, a numer- 
ical study of finite size systems is performed which is in exact agreement with the 
theoretical considerations. The general spectral properties of the non-abelian phase 
are considered for various finite toroidal systems. 
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1 Introduction 



Topological quantum computation [TPfSP] is certainly among the most ex- 
otic proposals for performing fault-tolerant quantum information processing. 
This approach has attracted considerable interest, since it is closely related to 
the problem of classifying topologically ordered phases in various condensed 
matter systems. The connection is provided by anyonic quasiparticles, which 
appear as states of topologically ordered systems with non-trivial statistical 
properties. Some of these anyon models can support universal quantum com- 
putation. Up to now, no complete classification of topological phases exists in 
terms of their physical properties or their computational power. This is due 
to the small number of analytically treatable models that exhibit topological 
behavior. The most studied arena is the celebrated fractional Quantum Hall 
effect appearing in a two dimensional electron gas when it is subject to 
a perpendicular magnetic field. 

Recently various two dimensional lattice models exhibiting topological behav- 
ior have been proposed [7f5|51[TU]|lllll2|[TB| that enjoy analytic tractability. One 
such lattice proposal is the honeycomb model introduced by Alexei Kitaev |lUj . 
It consists of a two dimensional honeycomb lattice with spins at its vertices 
subject to highly anisotropic spin-spin interactions. This model has several re- 
markable features. It is exactly solvable and can thus be studied analytically. 
For particular values of the couplings, the model can be mapped to Z2 gauge 
theory on a square lattice (the toric code), which supports abelian anyons. This 
anyon model has been employed for performing various quantum information 
tasks [2]. When one adds an external magnetic field, the model supports non- 
abelian Ising anyons. Even though neither model supports universal quantum 
computation, particular variations of the latter have been considered for this 
purpose [TlffT^] . One expects that when the couphngs of the honeycomb lat- 
tice model are varied, the system will undergo a phase transition between the 
abelian and non-abelian phases. The existence of the different phases is only 
argued in the original work [lOJ based on mathematical considerations and no 
rigorous presentation of the transition is provided. 

So far, the studies on Kitaev's honeycomb lattice model have concentrated 
on the abelian phase p^l7|18] . Here we present an extensive study of its 
spectral properties in the presence of an external magnetic field. Solving the 
model for various sparse vortex configurations gives us qualitative and quanti- 
tative results for the behavior of the spectrum in the non-abelian phase. The 
study includes the explicit demonstration of zero modes in the presence of well 
separated vortices and the lifting of the degeneracy due to their short-range 
interaction. These properties are subsequently connected to the properties of 
the Ising anyon model giving direct evidence that the low energy behavior of 
the non-abelian phase is indeed captured by this model. In addition, we con- 
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sider the stability of the different phases, which is of importance when one is 
interested in physically realizing the model [19] . The analytic calculations are 
supported by exact numeric diagonalizations of finite size systems. The ex- 
act agreement between the analytic and numeric solutions for these finite size 
systems is demonstrated and the effect of an external effective magnetic field 
on finite size systems is discussed. Our work generalizes the analysis in [T7] 
performed by one of the authors, where only the abelian phases in the limiting 
vortex-free and full-vortex cases were considered. 

The paper is organized as follows. In Section [2] we give an overview of the hon- 
eycomb lattice model. There we outline the analytic approach for solving the 
model for various vortex configurations by employing Majorana fermioniza- 
tion. Section [3] provides explicitly the analytic solution for the limiting cases of 
vortex-free and full-vortex configurations. These calculations are subsequently 
generalized to sparse vortex configurations. Section H] forms the main body of 
our work. There we analyze in detail the behavior of the spectrum in different 
parts of the phase space and study how it is modified due to the presence of 
vortices. A connection to the Ising anyon model is provided. In Section [5] we 
study exact numerical diagonalization of various finite size systems and show 
the equivalence with the analytic results. Final remarks and conclusions are 
given in Section O 



2 The spectrum of arbitrary vortex configurations 

2.1 The honeycomb lattice model 

We briefly review here the honeycomb lattice model and its analytical treat- 
ment as described in [10] . The model is defined on a honeycomb lattice A with 
spins residing at each site. The sites are bi-colored black and white such that 
A = As U A\Y, where A^ and A^y are two triangular sublattices. We shall 
consider the following Hamiltonian 

X— links y— links z— links 

where Jj,, Jy and are positive coupling strengths along the x-, y- and z- 
links, respectively, as shown in Figure [T](a). The three-spin interaction, or the 
i^-term, on the right hand side can be obtained from a perturbative expansion 
when we apply a weak (Zeeman) magnetic field of the form = h • cr. In this 
case K is given by ~ ^ g^^^ ^-j^jg model is assumed to approximate the 

one with a Zeeman term when K -C Jx, Jy, Jz- Only this third order term in 
the perturbative expansion will be of interest to us, since it is the lowest order 
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Fig. 1. The bi-colorable honeycomb lattice, (a) Depending on their orientation the 
links are labelled as x, y and z. (b) A single plaquette p with its sites enumerated, (c) 
Summation convention for each elementary unit cell. Solid arrows indicate nearest 
neighbor interactions along x- ,y- or z-links, whereas the dashed arrows indicate 
next to nearest neighbor interactions originating from the K-term [lOj. (d) The 
elementary unit cell with lattice basis vectors n^; and Uy. 

term breaking time reversal invariance [10]. The summations in the effective 
magnetic field term run over site triples such that every plaquette p contributes 
the six terms 



The sites of a single plaquette have been enumerated as shown in Figure [D^b). 
The Hamiltonian ([1]) commutes with the plaquette operators defined by 



w^ = alalalalalal, [1^;^ = /, (3) 

V 

where the product is taken over all plaquettes of a compact lattice and / 
is the identity operator. The eigenvalue ifp = — 1 is interpreted as having a 
vortex on plaquette p. The constraint in ([3]) implies that the vortices always 
come in pairs. Since are conserved quantities, one can fix the underlying 
vortex configuration and consider the Hamiltonian over this sector. This would 
not be possible when the usual Zeeman term were employed, since it does 
not commute with the plaquette operators. Thus, the magnetic field induces 
hopping of the vortices between neighboring plaquettes and their number is 
not necessarily conserved. 

The Hamiltonian can be diagonalized by representing the spin operators with 
Majorana fermions. Following [TOfTT] we introduce two fermionic modes a\ 
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and 02 residing at each lattice site. The corresponding Majorana fermions are 
given by 



We encode the spin at each site at the subspace where both of the fermionic 
modes are either empty or full. In terms of the four Majorana fermions, this 
means that we need to project down to a two dimensional subspace, the phys- 
ical space £, by employing the projector D = b^b^b^c, i.e. G £ -vv- = 
The representation of the spin matrices at site i is then given in terms 
of the Majorana fermions by = ib'^Ci, which satisfy the Pauli algebra when 
restricted in C (note that [Di, a"] = 0). It follows that 

crfcr" = -iiiijCiCj and cr^a^al = -iUikUjkCiCj, (5) 
where we have defined the operators 

Uij = ib^b"^, (uij = -Uji, uij = 1, ujj = Uij'j , (6) 

with a = x,y,z depending whether i E Kb and j G Kw are connected by a 
y- or 2;-link, respectively. Consequently, Hamiltonian ([1]) takes the form 

H = — AijCiCj, Aij = IJijUij + IK UjkUjk- (7) 

i jGA k 

The explicit appearance of the constraint D in ([5]) has been omitted, as we 
consider only operations in the physical subspace C. The couplings are given 
by Jij = J XI Jy or Jz- The summations in the Hamiltonian ([7]) are most con- 
veniently expressed pictorially (see Figure [U^c)). The solid lines correspond to 
nearest neighbor (the first term of Aij in ([7])) and the dashed lines to next to 
nearest neighbor summation (the second term of Aij in (^^). The antisymme- 
try of the m's in ([6]) is taken into account by using a convention such that one 
assigns an overall -|- (— ) to every term involving sites 2 G and j G 
when the arrow points from i to j {j to i). If two sites are not connected by 
an arrow the corresponding Aij element is zero. 

The plaquette operators ([3]) can be written in terms of the -u's as 

= n ^ ^ ^B, j e Aw, (8) 
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Fig. 2. The square lattice representation of the honeycomb lattice. Every vertex 
contains one black and one white site connected by a single z-link of the initial 
lattice. The basis vectors n^. and Uy (Figure [TJd)) become orthogonal unit vectors 
along the x- and y-links. 

where dp denotes the boundary of plaquette p. Also, one can check that 
[H, Uij] = 0. These observations imply that after performing the fermion- 
ization, the underlying vortex configuration can be fixed by specifying the 
eigenvalues Uij = ±1 of the operators Uij on every link of the model. The 
eigenvalue Uij = — 1 means that there is a string passing through the link 
ij that either connects two vortices or belongs to a loop. The locations p of 
these vortices are determined by the eigenvalues Wp = —1 of the plaquette 
operators. 

2.2 Solution for periodic vortex configurations 

Let us now consider in more detail the form of Hamiltonian ([7]) for various 
periodic vortex configurations and its diagonalization by using Fourier trans- 
form. Without affecting the physics of our system we shall deform the original 
honeycomb lattice to a square lattice by taking the length of 2;-links to zero 
and choosing the lattice basis vectors to be n^^ = (1,0) and n^^ = (0, 1). The 
resulting square lattice is shown in Figure [2J 

First we determine the unit cell of our periodic vortex lattice. The simplest 
possible choice contains a 2;-link of the honeycomb lattice, or in other words 
a single site on the square lattice. We refer to this choice of unit cell as the 
elementary cell. In order to employ Fourier transform in diagonalizing the 
Hamiltonian the underlying vortex configuration must be periodic with respect 
to the choice of the unit cell. For the elementary cell there is only one such 
configuration - the vortex-free configuration [10] . The full- vortex configuration, 
i.e. a vortex on every plaquette, can be solved by taking two elementary cells 
in which u alternates its sign along one direction [17]. However, our aim is 
to go beyond these two limiting cases and consider arbitrary sparse vortex 
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configurations. To do this, we define a generalized (M, A^)-unit cell, which 
contains MN elementary cells. On the square lattice the basis vectors take 
the simple form 



V, = Mn, = (M, 0), Wy = Nuy = (0, N). (9) 

An arbitrary site i on the original honeycomb lattice can be labeled by the 
triplet i (r, k, A), where r is a vector indicating the location of the unit cell, 
the index pair k = (m, n), 1 < m < M, 1 < n < A^, specifies a particular ele- 
mentary cell inside the generalized unit cell and A = 1, 2 is an index specifying 
whether the site belongs to or A^i/. If C denotes the number of unit cells, 
there are altogether 2MNC sites on the lattice. 

Using this notation the Hamiltonian ([7j) can be written as 



. iM,N) 2 

^=7E E E AAM(v)cA.(r)c^Kr + v), (10) 

^ r,v k,l \,fj,=l 

where the vector v is summed over all linear combinations of the lattice basis 
vectors ([9]). Since [H,u] = 0, the operators u appearing in Axk^ii dZl) have 
been replaced with their eigenvalues u = ±1. The antisymmetry of the oper- 
ators u is included in the summation convention (Figure [T](c)). The properly 
normalized Fourier transformation of the operators cxk is given by 



CAfc(r) 



tt/M 

2 r dp, 



tt/N 

■X f dp. 



Substituting this into (fTOl) we obtain the canonical form 



H- 



n/M 

1 f dp, 



n/N ( 
X I 0,Py 



2 J 27T/M J 

-n/M -tt/N 




Auip) A2(P) 

^21 (P) ^22 (P) 



Cl(P) 
\C2ip) 



where c{{p) = (4(i ^^^(p), . . ■ , cl^M,N)iP))^ ^nd AxM are matrices with ele- 
ments [ylA^(p)]fc/ = J2v'^^xkfii{'v)e~'^^'^ . This Hamiltonian is a generalization 
of the one obtained in [10] with the exception that the single entries of the 
2x2 Hamiltonian are replaced here with MN x MN matrices. 

The off-diagonal blocks correspond to nearest-neighbor interactions. The non- 
vanishing elements of the Hamiltonian for arbitrary (M, A^)-unit cells are given 
by 
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c|Ai2C2 = 2i{ +Uk,k Jz 



-\-n, 7 7 p»'5("-l)P-v!/ J r , ) 

^U'k,k-ny 'J yd '^l,k'^2,k-ny ) i 



(13) 



4>l2iCi = 2i( --Ujk^fe 4,ifcCi,fe 

-^^.,fe+n. J,e-^^(-^)P-^ 4,.Ci,,+„, (14) 

where the addition in the indices k = (m, n) is understood (m mod M, n mod 
A^) and 5{x) = 1 for x = and 5(,x) = otherwise. The diagonal blocks 
correspond to next-to-nearest neighbor couplings and are given by 





^-iS{n-N)p-Vy 




^l,k^l,k+ny 


k^k Tlx'^f^y 


gi(5(m-l)p-Vx Q-iS(n- 


■iV)p-V;, 


^\,k^l,k-nx+ny 


k ^ k~\~71yX 


g— ii5(m— M)p-Vx 




'^l,k'^l,k+nx 


k , k~{'Tix '^y 


g-ii5(m— M)p-Vj, gj(5(n 


-l)p-Vy 


^l,k^l,k+nx-ny 


^^k,k-nx 


^i5{m-l)p-Vx 




Ij/c 1)A/ TT-x 


k yk 


gi5(n-l)p-Vy 







(15) 



4^2202 = 





^—iS{n—N)p-Vy 




4,A;C: 


k ^k '^cc~\~'^y 


gi(5(m-l)p-Vx Q-i5(n- 


N)p-vy 


4,fc<^: 


\„,k+nx 
~'^k,k+nx 


^-i5(m-M)p-Vx 




4,fc'^: 


_yk+nx 

k,k-\-Tlx f^y 


^-i5{m—M)p-Vx giSin 


-1)P-V!/ 


4,itC; 


r^k 

kjk Tlx 


g?5(m-l)p-Vx 




4,jfcC; 


^ k^k "^y 


^iS{n-l)p-Vy 




4,fcC: 



(16) 



where we have used the short-hand notation ul. ^ = UkjUj^i. Both An and 
^422 are Hermitian, which can be checked by taking first Hermitian conjugates 
and subsequently shifting the indices accordingly. Likewise, one can check the 
relations, 
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All and A22 = -A^^. (17) 
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that guarantee the Hermiticity of A. 

The expressions derived above give the most general expression for the Hamil- 
tonian of the honeycomb lattice model. To proceed with the diagonalization, 
one needs to specify the underlying vortex configuration, i.e. the values of u 
on each link. Since all bi-colorable Hamiltonians have a double spectrum [10|, 
we know that the diagonalization of f|T2l) will give the general form 




where bi are MN fermionic operators and ei(p) are MN functions to be de- 
termined. The latter correspond to the eigenvalues ±ej(p) of the matrix A{p). 
Their exact form has to be calculated separately for each choice of unit cell 
and vortex configuration. The ground state and the first excited state corre- 
sponding to a particular vortex configuration are given by 

MN 

\9s)=ll n ^^(P)|0), |l) = 6l(po)M, (19) 

i=l —Tr<Px,Py<TT 

where |0) is a state with no Majorana fermions and po is the momentum 
minimizing the lowest lying eigenvalue ei(p), i.e. minp|ei(p)| = ei(po). It 
follows that the corresponding total ground state energy, E, and the fermion 
gap. A, are given by 




A = min|ei(p)|. (21) 



3 Analytic results at the thermodynamic limit 

In this section we present analytic solutions to the two limiting vortex configu- 
rations: the vortex-free and full- vortex configurations. Furthermore, we outline 
how the generalized unit cells can be used to study configurations, where the 
separation between two vortices is varied. This will be later used to study the 
behavior of the relative ground state energies and fermion gaps as the function 
of the vortex separation. 
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3.1 The vortex-free configuration 



The vortex-free configuration is achieved by setting 



«M = l,Vfc,/. (22) 

This configuration is periodic with respect to each ^-hnk and thus we can 
choose a (1, l)-unit cell (see Figure [21(a)). The off-diagonal, (IT^ . and diago- 
nal, (IT^ . elements are then given by 



Ai2(p) =22 (J, + J.e^P-^^ + Jye'^-'y) = 2/(p), 

^ii(p) =4ir (sin[p ■ (v^ - \y)] + sin(p ■ v^^) - sin(p ■ v^.)) = 5'(p)- 

Inserting these together with (fTTj) into (fT2|) . we obtain a 2 x 2 Hamiltonian 
which is diagonalized by introducing the fermionic operator 



6(p) = A ( C2(p) + r"^ ^i(p) 1 ' - '^^^^'^ 



/(p) ^^^7' (6(p)+^7(p))2 + |/(p)P 



where 



e(p) = ^/l/(p)P + ^?(p)^ 

|/(p)|^ = 4(J^ + + J^^ + 2 (J^. cosp^. + J^. Jj^ cos(p^. - py) + Jj^J^ cospj^)), 
g{Y)f = IGK"^ {sinpy - sinp^ + sm{p^ - py)f . 

The eigenvalues of the Hamiltonian are given by ±e(p), and thus the total 
ground state energy ( 120|) and the fermion gap (l2Ti) are given by 



} dp, r dpy e(p) 

^' = -J^J^^' ^^^^ 
Ao = mm|e(p)|. (24) 

These results agree with the ones obtained in [10] and [T7] . 
c?.^ The full-vortex configuration 



The full-vortex configuration can be obtained by choosing a (2, l)-unit cell 
and setting 
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Fig. 3. (a) The vortex-free configuration is created by setting u = 1 on all links. The 
configuration is periodic with respect to a (1, l)-unit cell, (b) The full-vortex con- 
figuration is created by alternating the value of the u's on y-links in the x-direction. 
We take ti = 1 on all x- and z-links. The configuration is periodic with respect to a 
(2, l)-unit cell. The solid squares denote the location of the vortices and the dashed 
lines indicate the strings along which u = —1. 



'1, 1) and I = k — n. 



otherwise. 



(25) 



Figure [3](b) illustrates the choice of unit cell for this case. Equations f[T3l) and 
( fTSi) become 



Ai2 = 2i 



and 



All = 2iK 



Jp-{Vy-Vx) 



^-tp vx -)- 1 -|- Q-^P-^y 



,ip vx — I — g^P-Vj, _ gip-(vi-Vj,) 
_^ip vy _j_ g-«p-Vy 



Inserting these into (fT2l) and diagonalizing the resulting 4x4 Hamiltonian we 
obtain the eigenvalues 



e(p) = ±2V/(p)±2V^?(p) 
where 



(26) 



/(p) = J2 + J2 ^ j2 ^ /^x\sm\p^ - py) + sin^ py + cos^Px), 

9{P) = JxJy COS^(Px ~Py) + J^Jz slu^ + JyJ^ COi? Py + 
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Jl siv? py + Jy cos^ + Jl svD?{px - Py) 

-{JxJy + JxJz + JyJz) sin(p^. - Py) sinpy cosp^. 



The analytic expressions of the corresponding eigenvectors are not presented 
here as they are too lengthy. When we set K = Q our results agree with the 
analytic calculations performed in [17] in the absence of the ii'-term. The total 
ground state energy fl20|) is now given by 



and the fermionic gap (|2T|) becomes 



A/-„ = min 
p 



(2J 



2^/(p)-2v/^7(p) 
c?..? Sparse vortex configurations 



We turn next to sparse vortex configurations in order to study the inter- 
actions between vortices. This is done by considering how the ground state 
energies and fermion gaps of 2-vortex configurations behave when the sepa- 
ration between vortices is varied. A configuration with two vortices separated 
by s plaquettes is created, for instance, by selecting an (M, A^)-unit cell and 
setting 



, —1, k = il < m < sA) and I = k — Uy 
Uk,i ={ y - - , ) y ^29) 

1, otherwise. 

Assuming M > N, vortex separations of s < M/2 can be studied. Figure H] 
illustrates the unit cell and the resulting vortex configuration. Since we are 
working on an infinite plane, ideally one would like to use a unit cell of infinite 
size in order to isolate the interaction between the vortices. However, the 
Hamiltonian f|T2|) grows polynomially in M and N and hence we restrict to 
considering a (20, 20)-unit cell. This choice allows separations s < 9, which is 
shown later to be sufficient to extract the expected asymptotic behavior when 
s ^ oo. The resulting Hamiltonians are sparse 800 x 800 matrices, which can 
not be treated analytically, but can be diagonalized numerically in a reasonable 
time using a tabletop computer. Using (!20l) and fl2T|) we can then calculate 
the total ground state energies and fermion gaps A|„ corresponding to 
2-vortex configurations with vortex separation s. 
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Fig. 4. (M, A^)-unit cell containing a pair of vortices separated by s plaquettes. 
s = means that the vortices occupy neighboring plaquettes. The solid squares in 
the plaquettes indicate the location of the vortices and the dashed line indicate the 
string along which u = —1 on all y-links. On all other links u = 1. 

4 Analysis of the spectrum 

The spectrum of Hamiltonian ([1]) can be characterized by two different types of 
energy gaps: fermionic gaps that characterize the energy levels of the spectrum 
above a fixed vortex configuration and vortex gaps that compare the ground 
state energies corresponding to different vortex configurations. We determine 
how the presence of vortices influences the fermionic gaps and, subsequently, 
how the phase space geometry is modified. We also determine the scaling of 
the ground state degeneracy, which is expected due to the presence of the 
Ising non-abelian vortices. Moreover, we carry out a study on the 2-vortex 
configuration energies as a function of the vortex separation, and subsequently 
determine the vortex gap, i.e. the energy required to excite a pair of free 
vortices. 



4-1 The fermion gap 

4-1.1 The phase space geometry in the presence of vortices 

First, we briefiy review the phase space of the vortex-free sector, which was 
studied in [TU]. It was shown that the honeycomb lattice model exhibits four 
distinct phases A^.Ay, and B for different values of the couplings such 
that the system is in the S-phase when all the inequalities \Jy\ + \Jz\ < \Jx\i 
\Jx\ + \Jz\ < \Jy \ and I Jx| + < \Jz\ are violated. The phase boundaries are 
given by the equalities and the phase A^ occurs when only | Jg| + \J^\ < \Ja\ 
holds and the other two inequalities are violated. The A^ phases are always 
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gapped for Jg, 7^ 0, K > and the vortices behave as Z2 x Z2 abehan 
anyons. On the other hand, the 5-phase is gapped only when K ^ and only 
there the vortices behave as non-abelian Ising anyons. The phase boundaries 
are the lines in the phase space where the fermion gap vanishes. Here we 
restrict to studying the transition, i.e. the behavior of the fermion gap between 
the A = Az (abelian) and the B (non-abelian) phases. Figure [5] illustrates 
the general phase space geometry where we've taken Jx + Jy + Jz = 1- For 
convenience we normalize from now on the couplings such that = 1 and 



FiguresE](a) andE](b) show the vortex- free, and full- vortex, 0281] . fermionic 
gaps plotted as functions of J for different values of K. Let us first consider 
the K = case. In the vortex-free configuration the gap vanishes at J = 1/2, 
in agreement with [10]. In the full- vortex configuration the gap persists till 
J = 1/ \/2 in line with the results derived in |T7] . There the phase bound- 
ary equalities for the full-vortex configuration were derived to take the form 
I Jgp + \ J'-^ = \ Ja\'- When i^' 7^ 0, in both cases the fermion gaps reappear 
and settle at a constant value once the system moves well into the B phase. 
However, the dependence of the gap magnitude on K is clearly different for the 
vortex-free and full-vortex configurations with the gap being much smaller in 
the latter case. Also, Figure[6]^b) shows that &XK = 1/5 the transition between 
the phases is shifted away from J = l/\/2. This implies that the magnitude 
of K can also affect the phase space geometry. 



In the limiting cases we observe that the boundary between the two phases 
depends on the underlying vortex configuration. To study how the fermionic 
energy gap interpolates in between these two extreme cases, we consider the 
fermion gaps of various sparse vortex configurations on small {MN < 12) unit 
cells. In all the cases the phase boundary falls into the region 1/2 < J < 1/^2 
when K = For very sparse configurations with low vortex density such as 
the 2- vortex configuration created by ( l29l) . the boundary is located very close 
to J = 1/2, whereas for more homogeneously distributed configurations with 
larger vortex density it tends towards J = 1/v^- However, when K > 
the phase boundary is in general shifted to larger J's such that for some 
configurations the boundary is located in the area J > 1/a/2. We attribute 
the shifting of the phase boundary to short-range vortex-vortex interactions, 
which are enhanced when K is increased. It is interesting to note that since the 
vortex density affects the phase space geometry in the region 1/2 < J < I/V2 
and K > 0, it could, in principle, be used as a tunable parameter, that induces 
a phase transition between the abelian and non-abelian phase. 
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Fig. 5. An illustration of the phase diagram with the four distinct phases Ax, Ay, A;, 
and B when + Jy + Jz = ^- We restrict to studying the transition only between 
A = Az and B phases and for convenience employ an alternative normalization such 
that J = Jx = Jy and = 1. The dashed line indicates the < J < 1 part of the 
phase space along which we study the system. For all vortex configurations at small 
K the phase boundary between A and B phases falls into the shaded area in the 
inner triangle. The limiting phase boundaries are given by the vortex-free (J = 1/2) 
and full- vortex ( J = l/\/2) configurations. 




J J 



(b) 

Fig. 6. (a) The vortex- free, (j24p . and (b) the full- vortex, ()28p . configuration fermion 
gaps for different values of K. In the vortex- free case the gap vanishes at J = 1/2 
for all values of K. In the full- vortex case the gap vanishes at J = l/\/2 for K = 0, 
but shifts to smaller J as is increased. 
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4-2 The fermion gap in the presence of vortices 



It is intriguing to study the behavior of the fermionic gap in the presence of 
only two vortices as a function of their distance. For that we consider again the 
2-vortex configurations with varying vortex separation s, which are created 
by setting the m's as given by (!29|) . Figure [7](a) shows the behavior of the 
corresponding fermion gap as the function of s in the abehan (J = 1/3) 
and in the non-abehan (J = 1) phase for several values of K. The behavior 
in the different phases is radically different. In the abelian phase the fermion 
gap is in practice insensitive to both s and K. In stark contrast to the abelian 
phase, the fermion gap in the non-abelian phase decreases exponentially with 
s vanishing completely for s > 2{3 Indeed, in Figure [7|(b) we plot the gaps to 
the two first excited states. 



and observe that tends exponentially to zero as s increases with the value 
at s = 9 being of order 10~^. This means that in the presence of two well sepa- 
rated vortices the ground state of the non-abelian phase is twofold degenerate. 
Moreover, the gap to the second excited state, A2v,2, is found to be insensitive 
to s and persist to arbitrary separations. 

It is known that the honeycomb lattice model can be mapped to a p-wave 
superconductor where the fermions live on the z-links [T61I25] . The vortices 
in the superconductor are assumed to be non-abelian Ising anyons [6], which 
are also predicted to appear as vortices in the honeycomb lattice model [TD]. 
In the presence of 2n well separated vortices the ground state should be 2"- 
fold degenerate, but when vortices are brought together, their interactions are 
predicted to lift the degeneracy Pl20ll24j . Our demonstration of the twofold 
degenerate ground state in the presence of two vortices is in agreement with 
this prediction. To verify that the degeneracy scales as 2" for the honeycomb 
lattice model, we consider in addition 4- and 6-vortex pairwise configurations 
where the vortex pairs are located on equally spaced rows of the (20, 20)- 
unit cell and the separation s of the vortices from each pair is simultaneously 
varied. In the 4-vortex case the fermion gaps to the four first excited states 
are given by 



The oscillations of the energy gap as a function of the vortex separation appear 
to be Friedel oscillations. 





(30) 
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(a) (b) 

Fig. 7. The gap behavior for a 2-vortex configuration as a function of the vortex 
separation s. (a) The fermion gap at the abehan (J = 1/3) and non-abehan 
(J = 1) phases for several values of K. (b) The two lowest lying excited states (pOl) 
in the non-abehan phase ( J = 1 and K = 1/5). In the non-abehan phase the ground 
state is a twofold degenerate for weU separated vortices. The degeneracy is lifted 
when the two vortices are brought close to each other. 





(a) (b) 

Fig. 8. The fermion gaps for some of the first excited states at J = 1 and K = 1/5 
as a function of the vortex separation s. (a) The fermion gaps (j3ip to the four 
first excited sates of the 4-vortex configuration. The first and second excited states 
are degenerate, (b) The fermion gaps ([32]) of the eight first excited states for the 
6- vortex configuration. 1st, 2nd and 3rd as well as 4th, 5th and 6th excited states 
remain degenerate at small s. 
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|1) = ^i(Po)l^s), A4^, = miiipo |ei(p)|, 

|2) = bl{po)\gs), A4^,,2 = miiipo |e2(p)|, ^^^^ 

|3) = bl{po)bl{po)\gs), A^^^s = minp,, |ei(p) + e2(p)|, 

|4) = bl{po)\gs), A4^,4 = minpo |e3(p)|. 



Similarly, the gaps to the eight first excited states for the 6-vortex configura- 
tion are 
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(32) 



Figures El^a) and (b) depict the behavior of the 4- and 6-vortex configuration 
fermion gaps (l3Ti) and ( l32l) . respectively, at J = 1 and ii" = 1/5 as a function 
of s. The degeneracy of the ground state as s — oo is four for 4- vortex and 
eight for 6-vortex configurations. This is exactly the predicted scaling. We also 
observe that at s < 2 the interaction does not completely hft the degeneracies. 
The ground state becomes non-degenerate for small s, but some of the states 
form degenerate bands. The splitting in energy between the bands is homoge- 
nous in the sense that for a specific s it costs the same amount of energy to 
move between the shifted states. We also observe that for all the considered 
vortex configurations the first non-vanishing gap as s — >■ cxo is always two units 
of energy above the ground state. As can be seen from Figure [6]^a), this is the 
case also for the vortex- free fermion gap, Aq, at J = 1 and K = 1/5. The 
observation Aq = A2„,2 = A4^^4 = Ag^.s suggests that the energy to excite a 
fermionic mode is insensitive to the underlying vortex configuration and de- 
pends only on J and K. We will adopt Aq to denote the energy to create a free 
fermionic excitation as opposed to an excited state at small s due to lifting of 
the ground state degeneracy. 
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4-2.1 Fermionic spectrum and the Ising anyon model 

Our results concerning the spectrum can be interpreted in the context of Ising 
anyon model, which is assumed to describe the low energy behavior of both 
the honeycomb lattice model and the p-wave paired superconductors PITU] . 
This model is spanned by three types of quasiparticles or sectors: the vacuum 
1, a fermion ip and a non-abelian a. The non-trivial fusion rules are given by 

ip X ip = 1, ip X a = a, axa = l + ip. (33) 

In the context of p-wave paired superconductors, 1 can be understood as the 
ground state condensate of Cooper pairs, ip a.s a. Bogoliubov quasiparticle and 
0" as a vortex In the honeycomb lattice model we take analogously 1 to 
be the ground state, to be a fermion mode bi in the spectrum (fT8|) and a to 
be a vortex living on a plaquette. 

An established method to study p-wave superconductor vortices is in terms of 
massless Majorana modes 7^ localized inside the vortex cores PI21II23II24] . Two 
Majorana modes can be combined to a fermion mode Zi = (74 + i7j+i)/\/2, 
which is carried by a pair of vortices located at i and i + Whether this mode 
is occupied or unoccupied corresponds to the two possible fusion outcomes of 
the a vortices - unoccupied mode corresponds to fusing to vacuum 1 whereas 
occupied mode means that the fusion will yield a ip. Since the occupation of 
these modes does not increase the energy of the system, they are known as 
zero modes, which appear in the spectra of systems supporting Ising anyons. 
The existence of n zero modes in the spectrum implies 2"-fold degenerate 
ground state. However, when the vortices are brought close to each other, 
the degeneracy should be lifted in a way that allows the determination of the 
fusion outcome [21I61I20] . 

The observed degeneracy in the presence of well separated vortices in the 
honeycomb lattice model (see Figures [71 E^a) and[8](b)) can be explained in 
terms of the "zero energy modes" in the spectrum. These modes do not strictly 
speaking have zero energy, but correspond instead to modes with the same 
finite energy as the ground state fl20|) . When n of these zero modes are present, 
we expect the diagonalized Hamiltonian (ITS]) to be of the form 




Here Q!|(p) are the n smallest eigenvalues that vanishes as the distance between 
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the vortices goes to infinity, i.e. linis^oo minp |a|(p)| = 0. We allow a|(p) to be 
finite at small s to account for the lifting of the ground state degeneracy. This 
is exactly what we obtain in the presence of vortices. Figures [7](b), [H](a) and 
MJo) show that for all s every occupied zero mode contributes equally to the 
energy splitting, which suggests ai(p) = a(p),Vz. This is reasonable, because 
occupied zero modes are interpreted as two a's fusing to a tp, and thus every 
mode at every s should contribute an equal energy proportional to the energy 
of a ip. 

The sphtting of the degenerate ground state in short ranges into degenerate 
bands spanned by states with the number of occupied zero modes can be used 
to extract information about how the vortices fuse. In particular, it is possible 
to identify the states with different fusion channels. Consider for instance the 
4-vortex configuration consisting of two well separated pairs whose fermion 
gap behavior is shown in Figure [Hl^b). The fusion rules fl55]) give 

axaxaxa = l + l + ilj + ilj, 

which mean that the four vortices may fuse to both vacuum 1 and to i{j in 
two distinct ways. These altogether four distinct fusion channels correspond 
to the fourfold degenerate ground state at large s. At small s there are three 
bands of different energy. The ground state corresponds to fusing both pairs 
into vacuum (no occupied zero modes), 

(cT X cr) X (cr X (t) — > 1 X 1 = 1. 

On the other hand, the non-degenerate band of two occupied zero modes 
corresponds also to the vacuum sector, but now such that both pairs will 
separately fuse to a ^Z', 

{a X a) X {a X a) ^ i/j X ip = 1. 

Even though this state belongs to the vacuum sector, they differ in energy 
because the two pairs are well separated from each other. This contrasts with 
the twofold degenerate band, which contains the states corresponding to the 
two fusion channels 

(a X a) X {a X a) 1 X ip = i/j and {a x a) x {a x a) tp x 1 = 

Both belong to the ip sector, but there is no energy splitting indicating which 
pair will fuse to ip and which to 1. This is actually only a feature of our con- 
struction where the vortices of the n pairs are always equidistantly separated 
by s plaquettes. As shown above, it is then only possible to deduce the total 
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topological sector of all the vortices, and some information about the global fu- 
sion channel by distinguishing between these states. However, if vortices from 
only one pair were brought close to each other, while the others were kept well 
separated, the interaction induced gap would correspond to a splitting of only 
one zero mode and the first excited state would be non-degenerate. Studying 
the splitting of each mode separately allows unambiguous determination of 
the global fusion channel. 

We also comment on the observation that the degree of degeneracy obtained 
here is 2". Usually one talks about creating vortices from vacuum, which means 
that the global sector of the system is fixed to 1. The degeneracy corresponding 
to 2n vortices is then 2"'"^, because when all the vortices are fused, one must 
obtain again the vacuum. Our observation of 2"-fold degeneracy means that 
we do not create vortices out of vacuum by fixing the m's over the unit cell. 
The overall sector of 2n vortices may be either 1 oi ip, but we have no prior 
information about it. This is reflected also on the identification of the fermion 
mode hi with a single if) excitation. If the overall sector was fixed, the V^'s 
should always appear in pairs. 

4.3 The vortex gap and the interaction energy 

We have studied the fermionic spectrum above a fixed background vortex con- 
figuration. Here we consider the energy spectrum of the vortices by studying 
how the ground state energy depends on the number of vortices and their sep- 
aration. A flux phase conjecture proven by Lieb [22] states that in the absence 
of an external field the energy minimum for honeycomb lattice is achieved 
with a vortex-free configuration. Even though we have an external magnetic 
field in our model, we assume this still holds when K is small. Under this 
assumption we define the vortex gap asymptotically by 

AE2^ = lim {E^J^ - MNEo), M > N, (35) 

M,N^oo 

which gives the energy to create a pair of vortices and drag them infinitely 
far from each other. Here Eq is the vortex-free ground state evaluated on a 
single plaquette fl23l) and denotes the total ground state energy of a vor- 
tex configuration on a (M, A^)-unit cell containing a pair of vortices separated 
by s plaquettes. The definition is given for a pair of vortices due to the con- 
straint ([3]) on the plaquette operators that demands vortices to come in pairs. 
Including the interaction energy in the vortex gap definition means that (135!) 
provides an estimate of the stability of the topological phase. To be precise, 
if the temperature of the system is well below the vortex gap, T « AE2V, 
spontaneous creation of stray vortex pairs will be exponentially suppressed. 
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1.2 




(a) (b) 

Fig. 9. The energy gap AE^ = E2^ — 20'^Eq between the vortex-free and the 2- vortex 
configuration with separation s on a (20,20)-unit cell in the (a) abelian (J = 1/3) 
and in the (b) non-abelian (J = 1) phase. 

To study how this definition applies to systems with periodic structure, we 
consider again the (20, 20)-unit cell with a pair of vortices separated by s 
plaquettes (l29l) . For a particular s, the vortex gap takes the form AE^ = 
— 20'^ Eq, which is plotted in Figure [9l The abelian phase (Figure M^a)) 
shows a very weak attractive short-range interaction between vortices, which 
is agreement with the high-order perturbation theory study of the abelian 
phase [18]. In the non-abelian phase (Figure [9]|^b)) we observe also an attractive 
interaction. However, there the interaction is strong with the magnitude being 
sensitive to the value of K. When K = the vortex gap exhibits only low 
amplitude oscillatory behavior as a function of s void of interaction signature, 
but as K increases and the system enters the non-abelian phase, the oscillatory 
behavior is suppressed and an attractive short-range interaction emerge^. In 
both phases the i^'-term increases the vortex energy, and in the non-abelian 
phase it is necessary to switch on the interaction. Although the interaction is 
very weak in the abelian phase, it exists even when K = 0. We observe that 
the interaction becomes negligible when s > 2, which is in accordance with 
the fermion gap behavior, which was attributed also to the interactions (See 
Figures Etb), (la) anddb)). 

Our definition contrasts with the vortex gap definition in [17], where it was 
defined as the difference between the total ground state energies of the full- 
vortex (!27j) and vortex-free fl23!l configurations. That definition neglected the 
interaction energy between vortices and hence we regard our definition (!35l) to 
provide a more accurate estimate of the energy required to excite the system. 
It should be emphasized that even though we observe a strong attractive inter- 
action between vortices, the plaquette operators ([3]) are constants of motion 



^ Like in Figure H^a), the physicality of the oscillatory behaviour for small K is 
unclear to us. 
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Fig. 10. The low energy spectrum of the non-abehan phase ( J = 1 and K = 1/5) as 
a function of vortex separation s. The plotted energies are with respect to Eq, the 
ground state energy of the vortex-free sector. The sohd hues are the total ground 
state energies and the dashed lines some of the lowest lying excited states above 
vortex- free (circle), 2- (square) and 4-vortex (diamond) configurations. 

of the Hamiltonian ([1]) and thus the vortices close to each other are not pulled 
together and annihilated spontaneously. This is not the case if the i^'-term was 
replaced by a Zeeman term. Then the vortices could hop and be annihilated 
if they do not have sufficient energy to overcome the attractive interaction. 
Therefore, we regard our asymptotic definition of the vortex gap fl5^ to pro- 
vide a reaUstic way of estimating the stability of the topological phase also in 
the presence of an external magnetic field. 



4-4 The low energy spectrum of the non-abelian phase 



Combining the studies on both the fermion and vortex gaps allows us to 
outline the low-energy spectrum of the non-abelian phase. In Figure [10] we 
plot some of the first excited fermionic states above the vortex-free, the 2- 
vortex and 4-vortex configurations as a function of the vortex separation s. 
All the energies are depicted with respect to the total ground state energy Eq 
of the vortex-free sector. The vortex-free sector is trivially insensitive to s and 
the corresponding low energy spectrum is characterized by the fermion gap Aq 
alone. In terms of the Ising anyon model these fermionic levels were identified 
with ip excitations. The 2- and 4-vortex sector ground states are separated 
from the vortex- free ground state by AE'^ and 2AE^, respectively, when the 
two pairs are well separated. At large s the first excited state above both 
of these configurations is given by the fermion gap to excite a free ip. Since 
this energy is constant on all vortex-configurations, these levels are located at 
AE"^ + Aq and 2AE'^ -|- Aq. As Aq > AE^ the low energy behavior consists 
only of vortices. Since the energy gap Aq to excite a ip does not depend 



23 



on the underlying vortex configuration, this observation generahzes to any 
configuration where the vortices are kept well separated. 

This simple spectral behavior is lost when the vortices are near each other. 
There the energy levels are shifted and the degeneracies are partially lifted. 
The smallest separation we can consider is s = 0, which corresponds to the 
vortices occupying neighboring plaquettes. However, when the vortices were 
superposed, they would fuse to the vacuum or to a fermion according to the 
fusion rules (|33ll . and the spectrum would correspond to the purely fermionic 
spectrum of the vortex-free sector. This can be connected to the lifting of the 
degeneracies at small s due to different number of occupied zero modes, i.e. 
different amount of V^'s in the fusion channels. We observe that at s = the 
ground states of both 2- and 4-vortex sectors (no occupied zero modes) tend 
towards the vortex-free ground state, whereas all the states with occupied zero 
modes tend to higher energies which correspond to fermionic excitations. This 
is in agreement with the predictions of the Ising anyon model. 



5 Numerical experiments on finite toroidal lattices 

In this section we present results of a numerical study of different finite size 
configurations. While the physics of these compactified systems is often com- 
plicated due to finite size effects, they can be used to directly compare nu- 
merical and analytical data. After a brief introduction into the methodology 
used in the numerical calculations, we present a comparison of numerical and 
analytical data and we show that they are in exact agreement, validating the 
presented theory. We continue with a more general numerical examination of 
finite toroidal systems. These studies go beyond the currently available an- 
alytical results and illustrate the non-trivial dependence of the low-energy 
spectrum on the i^T-term. 

5. 1 Systems of interest 

Our numerical experiments focus on calculating the low-energy spectral prop- 
erties of finite-size honeycomb lattice systems with the Hamiltonian given by 
([1]). We use a variety of toroidal systems which differ by the total number 
of spins and by lattice compactification. The size of the system varies from 
= 8 spins, which constitute an elementary unit of the model that can be 
compactified on a torus, to N = 24 spins. Though small scale (A^ < 18) com- 
putations require modest computing resources and are conveniently carried 
out using high level languages (e.g. Matlab), the dimensionality of the lattice 
Hilbert space scales exponentially with the number of spins (e.g. for 24 spins. 
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Fig. 11. Toroidal lattices used in numerical study. The diamond lattices (a) and (b), 
are symmetric with respect to exchange of x, y and z links. The rectangular lattices, 
(c), (d) and (e) are symmetric with respect to x and y-links only. 

the dimensionality of the Hilbert space is 2^*^ ^ 1.6 x 10^) and thus requires 
optimized parallel processing which will be discussed further below. In most 
cases the complexity of the computation can be reduced by taking into account 
the intrinsic symmetries of the model. These symmetries will be discussed in 
detail elsewhere [26]. We use two physically inequivalent types of finite lattice 
compactifications on a torus which we call 'diamond', (see Figures [TT] (a) and 
(b)) and 'rectangular', (see Figures fTTl (c). (d) and (e)). It can be easily seen 
that some nontrivial closed loops constructed within one compactification type 
correspond to open strings in the other type of the same size (for example the 
cases (a) and (c) in Figure [TT!). 

5.2 Methodology 

Our numerical methodology consists of three main components: generation of 
the Hamiltonian ([1]) and the plaquette operators Wp ([3]), their diagonalization 
and the analysis of the obtained eigenstates \n) and the eigenvalues En . 

For spin systems with (A^ > 18) we distribute the matrices over a number 
of different processors using the PETSc library [27|28II29] . The matrix loading 
routine requires the user to calculate the position and value of the non-zero ele- 
ments of a particular row of the matrix under consideration. Using this method 
we can easily store matrices of dimension 2^^. For basic linear algebra routines 
we use the Linear Algebra Wrapper (LAW) library [30]. For the numerical 
diagonalization of these matrices that are distributed across multiple proces- 
sors we use the SLEPc library [31| which is built upon PETSc. The software 
contains a number of exact diagonalization routines including the ARPACK 
Arnold! library [32||33] and an optimized implementation of Krylov-Shur algo- 
rithm [M1I35] . Using Krylov-Shur, and with the matrix distributed across 64 
processing nodes, SLEPc returns the lowest 10 energy eigensolutions of the 
full 24-spin system in under one hour. 

The aim of the numerical analysis is to classify the energy eigenvectors \n) ac- 
cording to their vortex configuration. Since all plaquettes commute with the 
Hamiltonian all energy eigenvectors \n) must satisfy Wp = {n\wp\n) = ±1 for 
all N/2 plaquettes on a torus. However, the relation UpW = I in implies 
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that there are only N/2 — 1 independent quantum numbers, {wi, ....,wiy/2-i} , 
for each vortex configuration sector. Therefore, in order to reduce the Hilbert 
space to particular vortex configuration, one must only impose N/2 — 1 con- 
straints. Since the imposition of each constraint reduces the dimension of the 
Hilbert space by a factor of 2 we see that there are only 2^/^~^ unique vortex 
configuration sectors, each with a Hilbert space dimension of 2^/^"*"^. 

In what follows we are only concerned with classifying the sectors according to 
the total number of vortices. For that we define the vortex counting operator 




The number of vortices corresponding to the eigenstate \n) is given by the 
expectation value {n\v\n). 

5. 3 Comparison of numerical and analytical data 

In this section we compare results obtained from exact numerical diagonaliza- 
tion with the solutions obtained by the method of Majorana fermionization 
outlined in Section 12. 2[ For the purpose of this comparison we analyze the 
= 8 spin diamond system (see Figure [TTlfa)). which corresponds to the 
(2, 2)-unit cell when we employ the Majorana fermionization. Placing this 
unit cell on a torus implies that all the couplings in the Hamiltonian ([7]) are 
now between sites belonging to the unit cell. In terms of the elements (1131) - 
flTBl) of the matrix A\k^i this means setting v^. = = everywhere. 

The (2, 2)-unit cell contains 12 links on which one must specify the values of the 
m's. This implies that there are altogether 2^^ = 4096 distinct ways to create 
vortex configurations. However, due to finite size effects there is no a priori 
way to tell which configuration of the u's will correspond to the lowest ground 
state energy Eq. We carry out a systematic investigation by diagonalizing the 
resulting 8x8 Hamiltonian for all the u configurations and find that there 
are in total 512, 3072 and 512 ways to create 0-vortex, 2-vortex and 4-vortex 
configurations, respectively. The lowest energy is found to correspond to the 
4-vortex configuration and the first excited states to vortex-free configurations. 

Using direct numerical diagonalization we find a ground state, which corre- 
sponds also to a 4-vortex configuration. In Figure [12] we plot some of the lowest 
lying excited states for (a) if = and (b) J = 1 and K >Q. The values agree 
with the analytical data to 14 decimal places. In addition, we observe using 
both numerical and analytic approach the level crossing due to the i^-term 
between 0-vortex and 4-vortex sectors (see Figure [TWb)). as well as that at 
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Fig. 12. Some of the lowest lying energy eigenvalues for the 8-spin diamond system 
(see Figure Eta)), (a) K = 0, (h) J = I and K > 0. 

J = 1 the vortex-free ground state is threefold degenerate. 
5.4 Numerical study of the K-term 

In this section we numerically calculate the effect of the K-term on the energy 
spectrum of the non-abelian B phase for three different finite size lattices. 
Analysis of the spectrum in all cases shows that the K-term is capable of 
inducing level crossings between eigenvectors from the same vortex configu- 
ration sector. This observation means that the K-term is, in the language of 
Section 12.21 capable of opening and closing of fermionic gaps. 

We first consider the 16-spin rectangular lattice shown in Figure [TT](d). The 
system is symmetric only with respect to reflection in the vertical (z-link) 
axis. At J = 1, = the ground state is four times degenerate, containing 
single 0- and 8- vortex states and two 4-vortex states. The two 4- vortex states 
are related by a lattice translation. In Figure [13] (a) we plot how the addition 
of the K-term affects the low energy spectrum. We observe a lifting of the 
ground state degeneracy, but the 4-vortex states remain still degenerate. We 
observe also that the K-term can induce spectral crossings between states from 
the same vortex configuration as seen in the double crossing of the 0-vortex 
ground states between K ^ 0.22 and K ^ 0.32. 

In Figure [T3](b) we consider the spectrum of 18-spin diamond lattice shown 
in Figure [TlTb). It is symmetric with respect to exchange of all x,y or z links 
and it contains an odd number of plaquettes. At the exact centre of the B 
phase (J = 1) the ground state contains three degenerate states belonging 
to the 0-vortex sector. We observe that the K-term does not lift the ground 
state degeneracy. This is to be expected as the K-term is also symmetric with 
respect to exchange of x,y oi z links. However, in general the K-term does 
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Fig. 13. Some of the lowest lying energy eigenvalues E above the ground state Eq 
for 16, 18 and 24-spin lattices at J = 1 and K > 0. (a) Spectrum of the 16 spin 
rectangular lattice of Figure [TlTdb (b) Spectrum of the 18 spin diamond lattice of 
Figure [TlTb). (c) Spectrum of the 24 spin rectangular lattice shown in Figure [TlTe) . 

affect the relative energy levels of non-degenerate states of the same vortex 
sector. This can be seen as the level crossings of the excited states. 

Finally, we investigate the 24-spin rectangular lattice shown in Figure [TlT e). 
We see that the lowest three states are non-degenerate and belong to the 
0-vortex sector. Again we observe the non-trivial behaviour of the spectrum 
as the function of K. Figure [T^ c) shows that the as K increases the gap 
between ground state and first excited state closes and a level crossing occurs 
&t K ^ 0.11. 



6 Conclusions 



We have presented an extensive analysis of the spectrum in Kitaev's honey- 
comb lattice model [lOj with the focus on the properties of the non-abelian 
regime. Due to the exact solvability of the model we were able to analyti- 
cally determine, qualitatively as well as quantitatively, the spectral behavior 
in the presence of vortices in the thermodynamic limit. This behavior was 
subsequently identified with the qualitative predictions derived in the context 
of a j9-wave superconductor, a model to which the honeycomb lattice model 
is known to be equivalent [TBf25] . The validity of our results is supported by 
exact numerical diagonalizations of various finite size lattice Hamiltonians. 
There is an exact agreement with the results obtained through the analytical 
methods. This is a strong validation of the employed analytical techniques and 
of the conclusions drawn from them. 

To be precise, our study allowed us to directly compare the spectral behavior 
in the abelian and in the non-abelian phases and extract characteristics that 
are unique to the non-abelian phase. The crucial difference is the ground 
state degeneracy in the non-abelian phase in the presence of well separated 
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vortices. We explain this in terms of zero modes in the spectrum and provide a 
direct verification that the number of zero modes present in a system with 2n 
vortices is n pT|23ll24j . The resulting 2"- fold degeneracy of the ground state is 
in agreement with the non-abelian character of the Ising non-abelian anyons. 
Furthermore, we observe directly the lifting of the ground state degeneracy 
when the vortices are brought close to each other and explain the lifting in 
terms of the fusion rules of the Ising vortices. The energy splitting at short 
ranges could in principle be used to distinguish the different possible fusion 
channels without the need to employ an interference procedure. Also, the fact 
that the information about the fusion outcome is a non-local property of the 
vortex pair is explicitly demonstrated as the degeneracy present at large s. 

Moreover, we have demonstrated that the phase boundary between the abelian 
and non-abelian phase depends on the underlying vortex configuration and 
that vortices are interacting in both phases. This attractive interaction is 
strong in the non-abelian and weak in the abelian phase. Also, the energy gap 
to excite a pair of vortices is considerably larger in the non-abelian phase. 
Another characteristic of the non-abelian phase is that the fermionic gap to 
excite free fermions, Aq, does not depend on the underlying vortex configura- 
tion as long as the vortices are kept well separated. This means that all the 
sectors of the non-abelian phase are equally stable with respect to fermionic 
excitations. Based on this we defined the vortex gap, Ai?2t;, in an asymptotic 
fashion. This gap provides a stability criterion for a particular sector. The 
combination of all these observations allowed us to outline the low energy 
spectrum of the non-abelian phase for configurations where the vortices are 
both well separated and close to each other. We observe that at large sepa- 
rations the spectral behavior consists only of vortices, which is in agreement 
with the prediction that the low energy behavior should be fully captured by 
the Ising anyon model. Understanding in detail the behavior of the energy 
spectrum is of importance to the proposed implementations of this model in 
the laboratory [19] . Our observation of the phase boundary dependence on the 
vortex density for particular values of Jq, and K could be of interest to these 
proposals. In particular, if an increase in temperature is accompanied by an 
increase in the number of vortices, then it would induce a transition from the 
non-abelian to the abelian phase. 

Finally, we have presented a numerical study of the effect of the i^'-term on the 
spectra of various finite size systems. These systems are too small to observe 
behavior similar to the thermodynamic limit, but one important qualitative 
similarity exists. We observe that the i^-term is capable of inducing level 
crossings of states belonging to the same vortex sector. This is the finite size 
equivalent to the opening and closing of fermionic gaps. A more detailed study 
on these finite size effects will be presented elsewhere [22] • 
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